subroutine Ljk(j,xlist0,xlist1,k, res)
! 作者：何光辉,
! 邮箱：flamehe@163.com
    !xlist0 denotes the original nodes, to construct the interpolation basis functions, 
    !e.g. the expanded Chebyshev nodes
    !xlist1 denotes the point to be interpolated, gauss duadrature points
    !the jth shape function of nodes xlist0 is to be returned, and then set x = xlist1[[k]]
    implicit none
    real(16) :: u, v, x, x0, xlist0(:), xlist1(:), res
    integer(4) :: length,m,j,k
    
    length = ubound(xlist0,1)
    x0 = xlist0(j); x = xlist1(k)
    u = 1q0; v = 1q0
    do m = 1, j-1
        u = u*(x - xlist0(m))
        v = v*(x0 - xlist0(m))
    end do
    do m = j+1, length
        u = u*(x - xlist0(m))
        v = v*(x0- xlist0(m))
    end do
    res = u/v
end subroutine Ljk
    